% 绘制简单的晶体
% 你或许可以用atomsk和ovito等工具来创建并展示更复杂的晶体
% CC-BY 4

clc
clear

function h = SHPERE(x,y,z,r)
  [X,Y,Z] = sphere(20);
  h = surf(r*X+x,r*Y+y,r*Z+z);
  set(h,'EdgeColor','none')
end

type = 1; %晶体种类： 0:SC, 1: BCC, 2:FCC, 3: HCP, 4: GRAPHITE
% 常温下部分金属的晶体种类参考:
% BCC: Fe
% FCC: Al, Cu, Au
% HCP: Ti

Nx = 4; %大小
Ny = 4;
Nz = 4;
R = 1;

figure()
axis equal
hold on

for i = 0:Nx-1
  for j = 0:Ny-1
    for k = 0:Nz-1

      isEdge = 0;
      if i == Nx - 1 || j == Ny - 1 || k == Nz-1
        isEdge = 1;
      endif

      switch type
        case 0
          a = 2;
          SPHERE(a*i,a*j,a*k,R);
        case 1
          a = 4/sqrt(3);
          SPHERE(a*i,a*j,a*k,R);
          if ~isEdge
            SPHERE(a*i+1/2*a,a*j+1/2*a,a*k+1/2*a,R);
          end
        case 2
          a = 4/sqrt(2);
          SPHERE(a*i,a*j,a*k,R);

          if ~(j == Ny - 1 || k == Nz-1) %xedge
            SPHERE(a*i,a*j+1/2*a,a*k+1/2*a,R);
          end
          if ~(i == Nx - 1 || k == Nz-1) %yedge
            SPHERE(a*i+1/2*a,a*j,a*k+1/2*a,R);
          end
          if ~(i == Nx - 1 || j == Ny-1) %zedge
            SPHERE(a*i+1/2*a,a*j+1/2*a,a*k,R);
          end
        case 3
          ax = 2;
          ay = 2*sqrt(3);
          az = sqrt(6)/3*4;
          SPHERE(ax*i, ay*j,az*k,R);

          if j ~= Ny-1
            SPHERE(ax*i + 1, ay*j + sqrt(3),az*k,R);
          end

          if i == 0 && j ~= Ny-1
            SPHERE(ax*i - 1, ay*j + sqrt(3),az*k,R);
          end

          if j ~= Ny-1 && k ~= Nz - 1
            if i ~= Nx - 1
              SPHERE(ax*i + 1, ay*j + 1/sqrt(3),az*k+az/2,R);
            end
            SPHERE(ax*i, ay*j + sqrt(3) + 1/sqrt(3) ,az*k+az/2,R);
          end

        case 4
          ax = 6;
          ay = 2*sqrt(3);
          az = 5;

          dx = 0;
          dy = 0;
          if mod(k,2) == 1
            dx=1;
            dy=sqrt(3);
          end

          SPHERE(ax*i+2+dx, ay*j+dy,az*k,R);
          SPHERE(ax*i+4+dx, ay*j+dy,az*k,R);
          if j ~= Ny - 1
            SPHERE(ax*i+5+dx, ay*j+sqrt(3)+dy,az*k,R);
            SPHERE(ax*i+1+dx, ay*j+sqrt(3)+dy,az*k,R);
          end
      end
    end
  end
end
drawnow

